import cupy as cp
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pyplot import MultipleLocator
import scipy.signal as sg

# 作图设置
plt.rcParams['axes.unicode_minus'] = False
plt.xticks(fontproperties='Times New Roman', fontsize=14)
plt.yticks(fontproperties='Times New Roman', fontsize=14)
# plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
# plt.rcParams.update({'font.size': 16})


class DataAnalysis:
    # 通用参数
    Re = 6371e3  # 地球半径
    c = 3e8  # 光速
    G = 6.67e-11  # 引力常量
    M = 5.965e24  # 地球质量

    def __init__(self, data, antenna_parameter, mission_parameter, radar_parameter, sample_parameter):
        # 天线参数读入
        self.Nx = antenna_parameter['Nx']
        self.Ny = antenna_parameter['Ny']
        self.W = antenna_parameter['W']
        self.dxRatio = antenna_parameter['dxRatio']
        self.dyRatio = antenna_parameter['dyRatio']

        # 雷达参数读入
        self.H = radar_parameter['H']
        self.f0 = radar_parameter['f0']
        self.PRF = radar_parameter['PRF']
        self.Br = radar_parameter['Br']
        self.N = radar_parameter['N']

        # 任务参数读入
        self.alpha_c = mission_parameter['alphaC']
        self.theta_sq_c = mission_parameter['thetaSqC']
        self.alpha_scan_min = mission_parameter['alphaScanMin']
        self.alpha_scan_max = mission_parameter['alphaScanMax']
        self.alpha_c_r = np.deg2rad(self.alpha_c)
        self.theta_sq_c_r = np.deg2rad(self.theta_sq_c)
        self.alpha_scan_max_r = np.deg2rad(self.alpha_scan_max)
        self.alpha_scan_min_r = np.deg2rad(self.alpha_scan_min)

        # 数据
        self.S0 = self.e2m(data['S0'])
        self.Sa = self.e2m(data['Sa'])
        self.w_opt = data['Wopt']

        # 采样手动设置
        self.R_sample_rate = sample_parameter['RSampleRate']
        self.f_sample_rate = sample_parameter['fSampleRate']

        # 雷达参数计算
        self.PRF_r = self.PRF * self.N
        self.pr = self.c / 2 / self.Br
        self.lamda = self.c / self.f0  # 波长
        self.v = (self.G * self.M / (self.H + self.Re)) ** 0.5  # 雷达平台速度
        self.f_eta_c = 2 * self.v * np.sin(self.theta_sq_c_r) / self.lamda  # 多普勒中心频率

        # 天线参数
        self.dx = self.dxRatio * self.lamda  # 天线x单元间隔
        self.dy = self.dyRatio * self.lamda  # 天线y单元间隔

        # 最大扫描角度参数计算
        self.alpha_ant_min = self.alpha_c - 90
        self.alpha_ant_max = self.alpha_c + 90
        self.alpha_ant_min_r = np.deg2rad(self.alpha_ant_min)
        self.alpha_ant_max_r = np.deg2rad(self.alpha_ant_max)
        self.theta_earth_r = np.arcsin(self.Re / (self.Re + self.H))  # 地球最大扫描角度
        self.theta_earth = np.rad2deg(self.theta_earth_r)

        # 主信号区
        self.theta_min_r = self.cal_theta(self.alpha_scan_min_r, self.theta_sq_c_r)
        self.theta_max_r = self.cal_theta(self.alpha_scan_max_r, self.theta_sq_c_r)
        self.theta_min = np.rad2deg(self.theta_min_r)
        self.theta_max = np.rad2deg(self.theta_max_r)

        # 雷达姿态
        self.A = np.arctan(np.tan(self.alpha_c_r) / np.cos(self.theta_sq_c_r))
        self.theta_yaw = np.arctan(np.sin(self.alpha_c_r) * np.tan(self.theta_sq_c_r))
        self.theta_pit = np.arcsin(np.cos(self.alpha_c_r) * np.sin(self.theta_sq_c_r))

        # 多普勒
        self.theta_bw_r = 0.886 / self.dyRatio / self.Ny / 4  # 波束宽度
        self.theta_bw = np.rad2deg(self.theta_bw_r)
        self.Bd = 2 * self.v * np.cos(self.theta_sq_c_r) * self.theta_bw_r / self.lamda  # 多普勒带宽

        self.R_min = self.theta2r_eta(self.theta_min_r)[0]  # 最近斜距
        self.R_max = self.theta2r_eta(self.theta_max_r)[0]  # 最远斜距

        # 计算模糊距离
        self.R0 = np.linspace(self.R_min, self.R_max, self.R_sample_rate)
        self.theta_nadir0_r = self.r2eta_theta(self.R0)[1]
        self.theta_nadir0 = np.rad2deg(self.theta_nadir0_r)

        # 模糊距离计算
        self.R_amb = []
        delta_r = self.c / 2 / self.PRF
        r_max = np.sqrt((self.Re + self.H) ** 2 - self.Re ** 2)

        # 测绘带到星下点（范围， 已包含镜像区）
        r = self.R0.copy()
        while True:
            r -= delta_r
            if r[0] < self.H and r[-1] < self.H:
                break
            elif r[0] < self.H < r[-1]:
                print('Warning: Nadir Echo')
                r[np.where(r < self.H)] = np.NAN
                self.R_amb.append(r.copy())
                break
            else:
                self.R_amb.append(r.copy())
        self.R_amb.reverse()

        # 测绘带到地球边界（已包含镜像区）
        r = self.R0.copy()
        while True:
            r += delta_r
            if r[-1] > r_max:
                break
            else:
                self.R_amb.append(r.copy())

    def r2eta_theta(self, r):  # 由斜距计算入射角，星下点角
        theta = np.arccos(np.round((r ** 2 + (self.Re + self.H) ** 2 - self.Re ** 2) / 2 / (self.Re + self.H) / r, 12))
        # alpha = np.arccos((r ** 2 + (Re + self.H) ** 2 - Re ** 2) / 2 / (Re + self.H) / r)
        eta = np.arcsin(np.round(np.sin(theta) * (self.Re + self.H) / self.Re, 12))
        # eta = np.arcsin(np.sin(alpha) * (Re + self.H) / Re)
        return eta, theta

    def cal_pattern_rotate(self, wmn):
        alpha = cp.linspace(self.alpha_ant_min, self.alpha_ant_max, 901)
        theta_sq = cp.linspace(-90, 90, 901)
        alpha_r = cp.deg2rad(alpha)
        theta_sq_r = cp.deg2rad(theta_sq)
        alpha, theta_sq = cp.meshgrid(alpha, theta_sq)
        alpha_r, theta_sq_r = cp.meshgrid(alpha_r, theta_sq_r)
        e = self.vector_2d(alpha_r, theta_sq_r)
        s = wmn.T @ e
        s = s.reshape(alpha.shape)
        s_db = abs(s)
        s_db = cp.asnumpy(s_db)
        s_db = s_db / np.nanmax(s_db)
        s_db = 20 * np.log10(s_db)
        s_db_raw = np.copy(s_db)
        s_db[np.where(s_db < -50)] = -50

        return s_db_raw, s_db
        # if if3d:
        #     ax = plt.axes(projection='3d')
        #     ax.plot_surface(alpha.get(), theta_sq.get(), s_db, rstride=1, cstride=1)
        #     ax.set_xlabel(r'$\alpha$')
        #     ax.set_ylabel(r'$\theta_{sq}$')
        #     ax.set_zlabel('S\\dB')
        #     ax.set_zlim([-50, 0])
        #     plt.figure()
        #
        # plt.figure()
        # plt.pcolor(alpha.get(), theta_sq.get(), s_db, cmap='jet', vmax=0, vmin=-50)
        # plt.colorbar()
        # plt.xlabel(r'$\alpha$')
        # plt.ylabel(r'$\theta_{sq}$')

    def cal_rar_a(self, ss, sa, wmn):
        wmn = wmn.reshape(self.Nx * self.Ny, 1)
        ss = np.sum(ss, axis=0)
        sa = np.sum(sa, axis=0)
        ss = wmn.T @ ss @ np.conj(wmn)
        sa = wmn.T @ sa @ np.conj(wmn)
        rar = (abs(sa) / abs(ss))
        rar_db = 10 * np.log10(rar)
        return rar_db

    def cal_rar_r(self, ss, sa, wmn):
        wmn = wmn.reshape(self.Nx * self.Ny, 1)
        ss = wmn.T @ ss @ np.conj(wmn)
        sa = wmn.T @ sa @ np.conj(wmn)
        rar = (abs(sa) / abs(ss)).reshape(self.R_sample_rate, )
        rar_db = 10 * np.log10(rar)
        # if if_plot:
        #     plt.plot(self.R0, rar_db)
        #     plt.xlabel('distance')
        #     plt.ylabel('RASR/dB')
        return rar_db

    def plot_asr_range(self):
        alpha_bg = np.linspace(-90, 90, 1081)
        theta_bg = np.linspace(-90, 90, 1081)
        alpha_bg, theta_bg = np.meshgrid(alpha_bg, theta_bg)
        background = np.ones(alpha_bg.shape) * 5
        plt.pcolor(alpha_bg, theta_bg, background, vmin=0, vmax=50, cmap='jet')

        # 多普勒频率
        f_eta_max = self.f_eta_c + self.Bd / 2
        f_eta_min = self.f_eta_c - self.Bd / 2
        f_eta0 = np.linspace(f_eta_min, f_eta_max, self.f_sample_rate)

        # 模糊距离
        r0 = self.R0.copy()
        r_amb = self.R_amb.copy()
        r_amb.append(r0.copy())
        r_amb = np.array(r_amb)

        for i in range(r_amb.shape[0]):
            for j in range(r_amb.shape[1]):
                eta, theta = self.r2eta_theta(r_amb[i, j])
                theta_sq = np.linspace(-theta, theta, int(2 * theta / np.pi * 180))
                alpha = self.cal_alpha(theta_sq, theta)
                plt.plot(np.rad2deg(alpha), np.rad2deg(theta_sq), color='y', linewidth=0.9)
                alpha = -alpha
                plt.plot(np.rad2deg(alpha), np.rad2deg(theta_sq), color='y', linewidth=0.9)
                print(f'{i}:{j}', end='\r')
        print('Range Ambiguity Finished')

        # 主信号区
        eta, theta = self.r2eta_theta(r0)
        theta_sq = self.doppler2theta_sq(f_eta0)
        for i in range(r0.shape[0]):
            alpha = self.cal_alpha(theta_sq, theta[i])
            plt.plot(np.rad2deg(alpha), np.rad2deg(theta_sq), color='r', linewidth=0.9)
        print('Main Signal Area Finish')
        plt.xlabel(r'$\alpha$', fontsize=14)
        plt.ylabel(r'$\theta_{sq}$', fontsize=14)

        num = 0
        # 测绘带及测绘带向上
        f_eta = f_eta0.copy()
        alpha = np.linspace(-90, 90, 1801)
        while True:
            theta_sq = self.doppler2theta_sq(f_eta)
            f_eta_next = f_eta + self.PRF_r
            theta_sq_next = self.doppler2theta_sq(f_eta_next)
            if np.any(np.isnan(theta_sq)):
                break
            theta_sq_interval = np.linspace(np.rad2deg(theta_sq[-1]), np.rad2deg(theta_sq_next[0]), 10)
            for i in range(theta_sq_interval.shape[0]):
                plt.plot(alpha, np.ones(1801,) * theta_sq_interval[i], color='b', linewidth='0.8')
            print(f'{num}', end='\r')
            # 停止条件
            num += 1
            f_eta = f_eta_next

        # 测绘带向下
        f_eta = f_eta0.copy()
        while True:
            f_eta_next = f_eta - self.PRF_r
            theta_sq = self.doppler2theta_sq(f_eta)
            theta_sq_next = self.doppler2theta_sq(f_eta_next)
            theta_sq_interval = np.linspace(np.rad2deg(theta_sq[0]), np.rad2deg(theta_sq_next[-1]), 10)
            for i in range(theta_sq_interval.shape[0]):
                plt.plot(alpha, np.ones(1801, ) * theta_sq_interval[i], color='b', linewidth='0.8')
            print(f'{num}', end='\r')
            f_eta = f_eta_next
            # 停止条件
            if np.any(np.isnan(theta_sq)):
                break
            num += 1

    def vector_2d(self, alpha, theta_sq):
        sx, sy = self.steer_2d(alpha, theta_sq)
        v = self.vector(sx, sy)
        shape = list(alpha.shape)
        shape.extend([self.Nx * self.Ny, 1])
        v = v.reshape(shape)
        return v

    def steer_2d(self, alpha, theta_sq):
        def phase(al, th):
            phx = ((cp.cos(self.A) * cp.cos(self.theta_yaw) +
                    cp.sin(self.A) * cp.sin(self.theta_yaw) * cp.sin(self.theta_pit)) * cp.sin(al) * cp.cos(th) +
                   (cp.cos(self.A) * cp.sin(self.theta_yaw) -
                    cp.sin(self.A) * cp.cos(self.theta_yaw) * cp.sin(self.theta_pit)) * cp.sin(th) -
                   cp.sin(self.A) * cp.cos(self.theta_pit) * cp.cos(al) * cp.cos(th)) * \
                  self.dx * 2j * cp.pi / self.lamda

            phy = (-cp.cos(self.theta_pit) * cp.sin(self.theta_yaw) * cp.sin(al) * cp.cos(th) +
                   cp.cos(self.theta_pit) * cp.cos(self.theta_yaw) * cp.sin(th) -
                   cp.sin(self.theta_pit) * cp.cos(al) * cp.cos(th)) * self.dy * 2j * cp.pi / self.lamda

            phz = ((cp.sin(self.A) * cp.cos(self.theta_yaw) -
                    cp.cos(self.A) * cp.sin(self.theta_yaw) * cp.sin(self.theta_pit)) * cp.sin(al) * cp.cos(th) +
                   (cp.sin(self.A) * cp.sin(self.theta_yaw) +
                    cp.cos(self.A) * cp.cos(self.theta_yaw) * cp.sin(self.theta_pit)) * cp.sin(th) +
                   cp.cos(self.A) * cp.cos(self.theta_pit) * cp.cos(al) * cp.cos(th))
            return phx, phy, phz

        phase_x, phase_y, phase_z = phase(alpha, theta_sq)

        shape_x = list(phase_x.shape)
        shape_y = list(phase_y.shape)
        shape_x.extend([self.Nx, 1])
        shape_y.extend([self.Ny, 1])
        sx = cp.zeros(tuple(shape_x), dtype='complex')
        sy = cp.zeros(tuple(shape_y), dtype='complex')
        p = cp.zeros(phase_x.shape, dtype='complex')
        for i in range(4):
            p += cp.exp(i * phase_y)
        for i in range(self.Nx):
            sx[..., i, 0] = cp.exp(i * phase_x)
        for i in range(self.Ny):
            sy[..., i, 0] = p * cp.exp(4 * i * phase_y)

        return sx, sy

    def e2m(self, s):
        m = np.zeros([s.shape[0], self.Nx * self.Ny, self.Nx * self.Ny], dtype='complex')
        t = s.copy()
        for i in range(self.Nx):
            m[:, self.Ny * i: self.Ny * (i + 1), :] = t
            if i < self.Nx - 1:
                t[:, :, self.Ny:] = t[:, :, 0: (self.Nx - 1) * self.Ny]
                t[:, :, 0: self.Ny] = np.transpose(np.conj(s[:, :, self.Ny * (i + 1): self.Ny * (i + 2)]), (0, 2, 1))
        return m

    def theta2r_eta(self, theta):  # 由星下点角计算斜距、入射角
        eta = np.arcsin(np.round(np.sin(theta) * (self.Re + self.H) / self.Re, 12))
        beta = eta - theta
        r = np.sin(beta) / np.sin(theta) * self.Re
        if type(r) == np.float64 and r == 0:
            r = self.H
        elif type(r) == np.array:
            r[np.where(beta == 0)] = self.H
        return r, eta

    def doppler2theta_sq(self, fd):
        return np.arcsin(self.lamda * fd / 2 / self.v)

    @staticmethod
    def cal_theta(alpha, theta_sq):  # input:radian output:radian
        return np.arccos(np.cos(alpha) * np.cos(theta_sq))

    @staticmethod
    def vector(sx, sy):
        sy = cp.transpose(sy, (0, 1, 3, 2))
        v = sx @ sy
        return v

    @staticmethod
    def cal_alpha(theta_sq, theta):
        return np.arccos(np.cos(theta) / np.cos(theta_sq))

    @staticmethod
    def cal_alpha_all(theta_sq, theta):
        theta_sq, theta = np.meshgrid(theta_sq, theta)
        return np.arccos(np.cos(theta) / np.cos(theta_sq)), theta_sq, theta

    @staticmethod
    def islr(p):
        # loc = np.where(p >= 10**(-0.3))
        # left = np.min(loc[0])
        # right = np.max(loc[0])
        # up = np.min(loc[1])
        # down = np.max(loc[1])
        # width = [right-left, down-up]
        # center = [int((right+left)/2), int((down+up)/2)]
        # locx = np.arange(center[0]-width[0], center[0]+width[0]+1,1)
        # locy = np.arange(center[1]-width[1], center[0]+width[1]+1,1)
        # locx, locy = np.meshgrid(locx,locy)
        p = p ** 2
        zero_x = sg.argrelmin(p, axis=0)
        zero_y = sg.argrelmin(p, axis=1)
        zeros = np.zeros(p.shape)
        zeros[zero_x[0], zero_x[1]] = 1
        zeros[zero_y[0], zero_y[1]] = 1
        center = np.where(p == 1)
        zero_x = np.where(zeros[:, center[1]] == 1)
        zero_y = np.where(zeros[center[0], :] == 1)
        zero_x = np.array(zero_x)
        zero_y = np.array(zero_y)
        width_x = np.min(abs(zero_x - center[0]))
        width_y = np.min(abs(zero_y - center[1]))
        locx = np.arange(center[0]-width_x, center[0]+width_x+1, 1)
        locy = np.arange(center[1]-width_y, center[1]+width_y+1, 1)
        locx, locy = np.meshgrid(locx,locy)
        return 10 * np.log10((np.sum(p) - np.sum(p[locx, locy])) / np.sum(p[locx, locy]))

    @staticmethod
    def pslr(p):
        px = np.max(p, axis=0)
        px = px[sg.argrelmax(px)]
        px.sort()
        py = np.max(p, axis=1)
        py = py[sg.argrelmax(py)]
        py.sort()
        return max(py[-2], px[-2])


if __name__ == '__main__':
    AntennaParameter = {
        'Nx': 22,
        'Ny': 16,
        'W': np.ones([22*16, 1]),
        'dxRatio': 0.65,
        'dyRatio': 0.64,
    }

    RadarParameter = {
        'H': 607e3,
        'f0': 1.26e9,
        'PRF': 3250,
        'Br': 80e6,
        'N': 1   # 通道数
    }

    MissionParameter = {
        'alphaC': 26,
        'thetaSqC': 15,
        'alphaScanMin': 25.06,
        'alphaScanMax': 27.33
    }

    # SampleParameter = {
    #     'RSampleRate': 22,
    #     'fSampleRate': 6
    # }

    SampleParameter = {
        'RSampleRate': 50,
        'fSampleRate': 5
    }
    Data = {
        'S0': np.load('S0_ARR.npy'),
        'Sa': np.load('Sa_ARR.npy'),
        'Wopt': np.load('W_ARR.npy')
    }

    data_analysis = DataAnalysis(Data, AntennaParameter, MissionParameter, RadarParameter, SampleParameter)

    Alpha = np.linspace(data_analysis.alpha_ant_min, data_analysis.alpha_ant_max, 361)
    Theta_sq = np.linspace(-90, 90, 361)
    Alpha_m, Theta_sq_m = np.meshgrid(Alpha, Theta_sq)

    pattern_raw, pattern = data_analysis.cal_pattern_rotate(cp.asarray(data_analysis.W))
    pattern_opt_raw, pattern_opt = data_analysis.cal_pattern_rotate(cp.asarray(data_analysis.w_opt))

    pattern_e = np.max(pattern_raw[100: 260, :], axis=0)
    pattern_opt_e = np.max(pattern_opt_raw[100: 260, :], axis=0)
    pattern_h = np.max(pattern_raw, axis=1)
    pattern_opt_h = np.max(pattern_opt_raw, axis=1)
    data_analysis.pslr(pattern_raw)

    rar_db = data_analysis.cal_rar_r(data_analysis.S0, data_analysis.Sa, data_analysis.W)
    rar_db_opt = data_analysis.cal_rar_r(data_analysis.S0, data_analysis.Sa, data_analysis.w_opt)

    rar_db_s = data_analysis.cal_rar_a(data_analysis.S0, data_analysis.Sa, data_analysis.W)
    rar_db_opt_s = data_analysis.cal_rar_a(data_analysis.S0, data_analysis.Sa, data_analysis.w_opt)
    print(f'{rar_db_s}dB {rar_db_opt_s}dB {rar_db_s-rar_db_opt_s}dB')

    # # fig0 = plt.figure(0)
    # # data_analysis.plot_asr_range()
    # # plt.savefig('./fig/fig0_300.tiff', format='tiff', dpi=300)
    # # plt.savefig('./fig/fig0_600.tiff', format='tiff', dpi=600)
    # # plt.show(block=True)
    #
    # fig1 = plt.figure(1)
    # plt.pcolor(Alpha_m, Theta_sq_m, pattern_raw, cmap='jet', vmax=0, vmin=-50)
    # plt.colorbar()
    # plt.xlabel(r'$\alpha$')
    # plt.ylabel(r'$\theta_{sq}$')
    # plt.savefig('./fig/fig1_300.tiff', format='tiff', dpi=300)
    # plt.savefig('./fig/fig1_600.tiff', format='tiff', dpi=600)
    #
    # fig2 = plt.figure(2)
    # ax = plt.axes(projection='3d')
    # ax.plot_surface(Alpha_m, Theta_sq_m, pattern_raw, rstride=1, cstride=1)
    # ax.set_xlabel(r'$\alpha$')
    # ax.set_ylabel(r'$\theta_{sq}$')
    # ax.set_zlabel(u'增益/dB')
    # plt.savefig('./fig/fig2_300.tiff', format='tiff', dpi=300)
    # plt.savefig('./fig/fig2_600.tiff', format='tiff', dpi=600)
    # # ax.set_zlim([-50, 0])
    #
    # fig3 = plt.figure(3)
    # plt.pcolor(Alpha_m, Theta_sq_m, pattern_opt_raw, cmap='jet', vmax=0, vmin=-50)
    # plt.colorbar()
    # plt.xlabel(r'$\alpha$')
    # plt.ylabel(r'$\theta_{sq}$')
    # plt.savefig('./fig/fig3_300.tiff', format='tiff', dpi=300)
    # plt.savefig('./fig/fig3_600.tiff', format='tiff', dpi=600)
    #
    # fig4 = plt.figure(4)
    # ax = plt.axes(projection='3d')
    # ax.plot_surface(Alpha_m, Theta_sq_m, pattern_opt_raw, rstride=1, cstride=1)
    # ax.set_xlabel(r'$\alpha$')
    # ax.set_ylabel(r'$\theta_{sq}$')
    # ax.set_zlabel(u'增益/dB')
    # plt.savefig('./fig/fig4_300.tiff', format='tiff', dpi=300)
    # plt.savefig('./fig/fig4_600.tiff', format='tiff', dpi=600)
    # # ax.set_zlim([-50, 0])
    #
    # fig5 = plt.figure(5)
    # plt.plot(Alpha, pattern_e, color='b', label=u'优化前')
    # plt.plot(Alpha, pattern_opt_e, color='r', label=u'优化后')
    # plt.legend()
    # plt.xlabel(r'$\alpha$')
    # plt.ylabel(u'增益/dB')
    # plt.savefig('./fig/fig5_300.tiff', format='tiff', dpi=300)
    # plt.savefig('./fig/fig5_600.tiff', format='tiff', dpi=600)
    #
    # fig6 = plt.figure(6)
    # plt.plot(Theta_sq, pattern_h, color='b', label=u'优化前')
    # plt.plot(Theta_sq, pattern_opt_h, color='r', label=u'优化后')
    # plt.legend()
    # plt.xlabel(r'$\theta_{sq}$')
    # plt.ylabel(u'增益/dB')
    # plt.savefig('./fig/fig6_300.tiff', format='tiff', dpi=300)
    # plt.savefig('./fig/fig6_600.tiff', format='tiff', dpi=600)
    #
    # fig7 = plt.figure(7)
    # plt.plot(data_analysis.R0 / 1000, rar_db, color='b', label=u'优化前')
    # plt.plot(data_analysis.R0 / 1000, rar_db_opt, color='r', label=u'优化后')
    # plt.legend()
    # plt.xlabel(u'斜距/km')
    # plt.ylabel('ASR/dB')
    # plt.savefig('./fig/fig7_300.tiff', format='tiff', dpi=300)
    # plt.savefig('./fig/fig7_600.tiff', format='tiff', dpi=600)
    #
    # fig8 = plt.figure(8)
    # w = data_analysis.w_opt
    # w = w / abs(np.max(w))
    # w = abs(w)
    # x = np.arange(data_analysis.Nx).reshape(-1, 1)
    # y = np.arange(data_analysis.Ny).reshape(-1, 1)
    # x, y = np.meshgrid(x, y)
    # w = w.reshape(x.shape)
    # x, y, w = x.ravel(), y.ravel(), w.T.ravel()
    # y_major_locator = MultipleLocator(4)
    # ax = plt.axes(projection='3d')
    # # ax.plot_surface(x, y, w, rstride=1, cstride=1, cmap='jet')
    # ax.bar3d(x, y, 0, 0.7, 0.7, w, cmap='jet')
    # ax.set_zlim([0, 1])
    # ax.set_xlabel(u'俯仰向单元')
    # ax.set_ylabel(u'方位向单元')
    # ax.yaxis.set_major_locator(y_major_locator)
    # ax.set_zlabel('|w|')
    # plt.savefig('./fig/fig8_300.tiff', format='tiff', dpi=300)
    # plt.savefig('./fig/fig8_600.tiff', format='tiff', dpi=600)
    #
    # fig9 = plt.figure(9)
    # w = data_analysis.w_opt
    # w = w / abs(np.max(w))
    # w = abs(w)
    # x = np.arange(data_analysis.Nx).reshape(-1, 1)
    # y = np.arange(data_analysis.Ny).reshape(-1, 1)
    # x, y = np.meshgrid(x, y)
    # w = w.reshape(x.shape)
    # y_major_locator = MultipleLocator(4)
    # ax = plt.axes(projection='3d')
    # ax.plot_surface(x, y, w, rstride=1, cstride=1, cmap='jet')
    # ax.set_zlim([0, 1])
    # ax.set_xlabel(u'俯仰向单元')
    # ax.set_ylabel(u'方位向单元')
    # ax.yaxis.set_major_locator(y_major_locator)
    # ax.set_zlabel('|w|')
    # plt.savefig('./fig/fig9_300.tiff', format='tiff', dpi=300)
    # plt.savefig('./fig/fig9_600.tiff', format='tiff', dpi=600)
    #
    # fig10 = plt.figure(10)
    # y_major_locator = MultipleLocator(2)
    # x_major_locator = MultipleLocator(2)
    # plt.pcolor(x, y, w, cmap='jet', vmax=1, vmin=0)
    # plt.colorbar()
    # plt.xlabel(u'俯仰向单元')
    # plt.ylabel(u'方位向单元')
    # ax = plt.gca()
    # ax.xaxis.set_major_locator(x_major_locator)
    # ax.yaxis.set_major_locator(y_major_locator)
    # plt.savefig('./fig/fig10_300.tiff', format='tiff', dpi=300)
    # plt.savefig('./fig/fig10_600.tiff', format='tiff', dpi=600)
    #
    # fig0 = plt.figure(0)
    # data_analysis.plot_asr_range()
    # plt.savefig('./fig/fig0_300.tiff', format='tiff', dpi=300)
    # plt.savefig('./fig/fig0_600.tiff', format='tiff', dpi=600)
    # plt.xticks(fontproperties='Times New Roman', fontsize=14)
    # plt.yticks(fontproperties='Times New Roman', fontsize=14)
    # plt.show(block=True)

    fig1 = plt.figure(1)
    plt.pcolor(Alpha_m, Theta_sq_m, pattern_raw, cmap='jet', vmax=0, vmin=-50)
    plt.colorbar()
    plt.xlabel(r'$\alpha$', fontsize=14)
    plt.ylabel(r'$\theta_{sq}$', fontsize=14)
    plt.xticks(fontproperties='Times New Roman', fontsize=14)
    plt.yticks(fontproperties='Times New Roman', fontsize=14)
    plt.savefig('./fig/fig1_300.jpg', format='jpg', dpi=300)
    plt.savefig('./fig/fig1_600.jpg', format='jpg', dpi=600)

    fig2 = plt.figure(2)
    ax = plt.axes(projection='3d')
    ax.plot_surface(Alpha_m, Theta_sq_m, pattern_raw, rstride=1, cstride=1, cmap='jet')
    ax.set_xlabel(r'$\alpha$', fontsize=14)
    ax.set_ylabel(r'$\theta_{sq}$', fontsize=14)
    ax.set_zlabel(u'增益/dB', fontsize=14)
    # plt.tick_params(labelsize=14)
    plt.tight_layout()
    plt.savefig('./fig/fig2_300.jpg', format='jpg', dpi=300)
    plt.savefig('./fig/fig2_600.jpg', format='jpg', dpi=600)
    # ax.set_zlim([-50, 0])

    fig3 = plt.figure(3)
    plt.pcolor(Alpha_m, Theta_sq_m, pattern_opt_raw, cmap='jet', vmax=0, vmin=-50)
    plt.colorbar()
    plt.xlabel(r'$\alpha$', fontsize=14)
    plt.ylabel(r'$\theta_{sq}$', fontsize=14)
    plt.xticks(fontproperties='Times New Roman', fontsize=14)
    plt.yticks(fontproperties='Times New Roman', fontsize=14)
    plt.savefig('./fig/fig3_300.jpg', format='jpg', dpi=300)
    plt.savefig('./fig/fig3_600.jpg', format='jpg', dpi=600)

    fig4 = plt.figure(4)
    ax = plt.axes(projection='3d')
    ax.plot_surface(Alpha_m, Theta_sq_m, pattern_opt_raw, rstride=1, cstride=1, cmap='jet')
    ax.set_xlabel(r'$\alpha$', fontsize=14)
    ax.set_ylabel(r'$\theta_{sq}$', fontsize=14)
    ax.set_zlabel(u'增益/dB', fontsize=14)
    # plt.tick_params(labelsize=14)
    plt.tight_layout()
    plt.savefig('./fig/fig4_300.jpg', format='jpg', dpi=300)
    plt.savefig('./fig/fig4_600.jpg', format='jpg', dpi=600)
    # ax.set_zlim([-50, 0])

    fig5 = plt.figure(5)
    plt.plot(Alpha, pattern_e, color='b', label=u'优化前')
    plt.plot(Alpha, pattern_opt_e, color='r', label=u'优化后')
    plt.legend(fontsize=14)
    plt.xticks(fontproperties='Times New Roman', fontsize=14)
    plt.yticks(fontproperties='Times New Roman', fontsize=14)
    plt.xlabel(r'$\alpha$', fontsize=14)
    plt.ylabel(u'增益/dB', fontsize=14)
    plt.savefig('./fig/fig5_300.jpg', format='jpg', dpi=300)
    plt.savefig('./fig/fig5_600.jpg', format='jpg', dpi=600)

    fig6 = plt.figure(6)
    plt.plot(Theta_sq, pattern_h, color='b', label=u'优化前')
    plt.plot(Theta_sq, pattern_opt_h, color='r', label=u'优化后')
    plt.legend(fontsize=14)
    plt.xticks(fontproperties='Times New Roman', fontsize=14)
    plt.yticks(fontproperties='Times New Roman', fontsize=14)
    plt.xlabel(r'$\theta_{sq}$', fontsize=14)
    plt.ylabel(u'增益/dB', fontsize=14)
    plt.savefig('./fig/fig6_300.jpg', format='jpg', dpi=300)
    plt.savefig('./fig/fig6_600.jpg', format='jpg', dpi=600)

    fig7 = plt.figure(7)
    plt.plot(data_analysis.R0 / 1000, rar_db, color='b', label=u'优化前')
    plt.plot(data_analysis.R0 / 1000, rar_db_opt, color='r', label=u'优化后')
    plt.legend(fontsize=14)
    plt.xticks(fontproperties='Times New Roman', fontsize=14)
    plt.yticks(fontproperties='Times New Roman', fontsize=14)
    plt.xlabel(u'斜距/km', fontsize=14)
    plt.ylabel('ASR/dB', fontsize=14)
    plt.savefig('./fig/fig7_300.jpg', format='jpg', dpi=300)
    plt.savefig('./fig/fig7_600.jpg', format='jpg', dpi=600)

    # fig8 = plt.figure(8)
    # w = data_analysis.w_opt
    # w = w / abs(np.max(w))
    # w = abs(w)
    # x = np.arange(data_analysis.Nx).reshape(-1, 1)
    # y = np.arange(data_analysis.Ny).reshape(-1, 1)
    # x, y = np.meshgrid(x, y)
    # w = w.reshape(x.shape)
    # x, y, w = x.ravel(), y.ravel(), w.T.ravel()
    # y_major_locator = MultipleLocator(4)
    # ax = plt.axes(projection='3d')
    # # ax.plot_surface(x, y, w, rstride=1, cstride=1, cmap='jet')
    # ax.bar3d(x, y, 0, 0.7, 0.7, w, cmap='jet')
    # ax.set_zlim([0, 1])
    # ax.set_xlabel(u'俯仰向单元')
    # ax.set_ylabel(u'方位向单元')
    # ax.yaxis.set_major_locator(y_major_locator)
    # ax.set_zlabel('|w|')
    # plt.savefig('./figM/fig8_300.jpg', format='jpg', dpi=300)
    # plt.savefig('./figM/fig8_600.jpg', format='jpg', dpi=600)

    # fig9 = plt.figure(9)
    w = data_analysis.w_opt
    w = w / abs(np.max(w))
    w = abs(w)
    x = np.arange(data_analysis.Nx).reshape(-1, 1)
    y = np.arange(data_analysis.Ny).reshape(-1, 1)
    x, y = np.meshgrid(x, y)
    w = w.reshape(x.shape)
    # y_major_locator = MultipleLocator(4)
    # ax = plt.axes(projection='3d')
    # ax.plot_surface(x, y, w, rstride=1, cstride=1, cmap='jet')
    # ax.set_zlim([0, 1])
    # ax.set_xlabel(u'俯仰向单元')
    # ax.set_ylabel(u'方位向单元')
    # ax.yaxis.set_major_locator(y_major_locator)
    # ax.set_zlabel('|w|')
    # plt.savefig('./figM/fig9_300.jpg', format='jpg', dpi=300)
    # plt.savefig('./figM/fig9_600.jpg', format='jpg', dpi=600)

    fig10 = plt.figure(10)
    y_major_locator = MultipleLocator(2)
    x_major_locator = MultipleLocator(2)
    plt.pcolor(x, y, w, cmap='jet', vmax=1, vmin=0)
    plt.colorbar()
    plt.xlabel(u'俯仰向单元', fontsize=14)
    plt.ylabel(u'方位向单元', fontsize=14)
    plt.xticks(fontproperties='Times New Roman', fontsize=14)
    plt.yticks(fontproperties='Times New Roman', fontsize=14)
    ax = plt.gca()
    ax.xaxis.set_major_locator(x_major_locator)
    ax.yaxis.set_major_locator(y_major_locator)
    plt.tight_layout()
    plt.savefig('./fig/fig10_300.jpg', format='jpg', dpi=300)
    plt.savefig('./fig/fig10_600.jpg', format='jpg', dpi=600)

    print('Program End')
    # print('Program End')